int toResultIdx(int row, int col)
{
	if(row<col)
	{
		return (col-1)*col/2+row;
	}
	else
	{
		return (row-1)*row/2+col;
	}
}

__kernel void test(
	__global float *data,   //D*N
	__global float *result, //(blocksize-1)*blocksize/2, for first stage
	__global unsigned int *row_counter,
	__global unsigned int *cross_counter,
	__local float *block1,
	__local float *block2,
	__local unsigned int *current_row,
	unsigned int N,
	unsigned int D,
	unsigned int blocksize)
	
{
	//int igrid = get_global_id(0);
	int lsize = get_local_size(0);
	int lid = get_local_id(0);
	int gid;
	float diag;
	float temp;
	float a;
	float b;
	int idx;
	int globidx;
	int i;
	int j;
	
	
	//for now, lsize must be >= than blocksize
	//preferably they are equal
	
	//load first block into local memory
	//this one is slightly slow because not coalesced (data is column-major, block1 is row-major)
	//but we only do it once
	/* column-major code
	for(int i=0; i<D*blocksize; i+=lsize)
	{
		if(lid+i<blocksize*D)
		{
			block1[i]=data[gid*D*blocksize+lid+i];
		}
	}
	*/
	if(lid==0)
	{
		gid=atom_inc(row_counter);
		*current_row=gid;
	}
	barrier(CLK_LOCAL_MEM_FENCE);
	gid=*current_row;
	barrier(CLK_LOCAL_MEM_FENCE);

	while(gid<N/blocksize)
	{
		globidx=gid*blocksize+lid;
		
	
	diag=0;
	temp=0;
	//row major load, and compute square terms at the same time
	
	for(i=0; i<D; i++)
	{
		if(lid<blocksize)
		{
			temp=data[blocksize*gid+N*i+lid];
			block1[i*blocksize+lid]=temp;
			diag+=temp*temp;
		}
	}
	
	barrier(CLK_LOCAL_MEM_FENCE);
	
	
	
	
	
	
	//main loop for computing cross terms
	for(int col=N/blocksize - 1; col>gid; col--)
	{
		if(col>gid)
		{
			//load another block into local mem
			
			//for(i=0; i<D*blocksize; i+=lsize)
			//{
			//	if(lid+i<blocksize*D)
			//	{
			//		block2[lid+i]=data[col*D*blocksize+lid+i];
			//	}
			//}
			
			
			//row-major
			for(int j=0; j<D; j++)
			{
				if(lid<blocksize)
				{
					block2[j*blocksize+lid]=data[blocksize*col+N*j+lid];
				}
			}
			
			barrier(CLK_LOCAL_MEM_FENCE);
	
			//compute the dot product of each column.
			
			if(lid<blocksize)
			{
				
				for(int k=0; k<blocksize; k++)
				{
					temp=0;
					//row-major, fast
					
					for(i=0; i<D; i++)
					{
						a=block2[i*blocksize+k];
						b=block1[i*blocksize+lid];
						temp+=a*b;
					}
					
					idx=toResultIdx(globidx,col*blocksize+k);
					result[idx]=diag-2*temp;
				}
				
				
				
			}
			
		if(lid==0)
		{
			atom_inc(&cross_counter[col]);
		}
		barrier(CLK_LOCAL_MEM_FENCE);
		
		}
		//barrier(CLK_LOCAL_MEM_FENCE);
		
		
		
	}
	
	
	
	//compute cross terms adjacent to the diagonal
	for(j=1; j<blocksize; j++)
	{
		temp=0;
		//row-major, fast
		for(int i=0; i<D; i++)
		{
			a=block1[i*blocksize+j];
			b=block1[i*blocksize+lid];
			temp+=a*b;
		}
	
		if(lid<j)
		{
			idx=toResultIdx(globidx,j+gid*blocksize);
			result[idx]=diag-2*temp;
		}
	}
	
	
	
	//wait for others to finish column above me
	 
	if(lid==0)
	{
		while(cross_counter[gid]<gid){}
	}
	
	barrier(CLK_LOCAL_MEM_FENCE);
	
	/*
	int offset=(globidx-1)*globidx/2;
	if(lid<blocksize)
	{
		for(int i=0; i<globidx; i++)
		{
			result[i+offset]=result[i+offset]+diag;
		}
	
	}
	*/
	
	int offset;
	
	block1[lid]=diag;
	barrier(CLK_LOCAL_MEM_FENCE);
	
	for(i=0; i<blocksize; i++)
	{
		diag=block1[i];
		//barrier(CLK_LOCAL_MEM_FENCE);
		idx=lid;
		globidx=gid*blocksize+i;
		offset=toResultIdx(0,globidx);
		
			for(int j=0; j<=gid; j++)
			{
				if(idx<globidx)
				{
					result[offset+idx]=result[offset+idx]+diag;
				}
				idx+=blocksize;
			}
		
	}
	
	
	if(lid==0)
	{
		gid=atom_inc(row_counter);
		*current_row=gid;
	}
	barrier(CLK_LOCAL_MEM_FENCE);
	gid=*current_row;
	//barrier(CLK_LOCAL_MEM_FENCE);

	
	}
	
	
	
}